function [GlobKe,GlobBe]=CalculatePlaneKeBe(NodePt,PlaneNd)
global charge

GlobKe=zeros(length(NodePt));
GlobBe=zeros(length(NodePt),1);

PlaneSize=size(PlaneNd);
for n=1:PlaneSize(1)
    Coord=zeros(3,2);
    for np=1:3
        Coord(np,1)=NodePt(PlaneNd(n,np+1),2);
        Coord(np,2)=NodePt(PlaneNd(n,np+1),3);
    end
     
    Ke=zeros(3,3);
    Be=zeros(3,1);
    GsNum=3;
    for inp=1:GsNum
        [zeta,weight]=GetTriGaussIntgInfo(GsNum,inp);
        [Nshape,dNdzeta]=GetTriNshape(3,zeta);
        Jm=dNdzeta*Coord;
        GradN=Jm\dNdzeta;
        Jdet=det(Jm);
        Ke=Ke+GradN'*GradN*weight*Jdet;
    end
    
    for ni=1:3;
        FlagNi=PlaneNd(n,ni+1);
        GlobBe(FlagNi)=GlobBe(FlagNi)+Be(ni);
        for nj=1:3
            FlagNj=PlaneNd(n,nj+1);
            GlobKe(FlagNi,FlagNj)=GlobKe(FlagNi,FlagNj)+Ke(ni,nj);
        end
    end
end